#############################################
## PROGRAM NAME: 203f_res_muni_noimp.R     ##
## AUTHOR: MATT MLECZKO                    ##
## INPUTS:                                 ##
##         201_af_muni_fin.Rda             ##
##         202_muni_results_supp_1.qs -    ##
##         202_muni_results_supp_100.qs    ##
##                                         ##
## OUTPUTS:                                ##
##         203f_m1_muni_noimp.Rda -        ##
##         203f_m9_muni_noimp.Rda          ##
##                                         ## 
## PURPOSE: Plot muni no imp results       ##
##                                         ##
## LIST OF UPDATES:                        ##
#############################################

#log <- file(# USER DEFINED PATH AND FILE NAME HERE #)
#sink(log, append=TRUE)
#sink(log, append=TRUE, type="message")

## load libraries ##
library(tidyverse) 
library(qs)
library(stringr)

data_path <- # USER DEFINED PATH HERE #
setwd(data_path)

load("201_af_muni_fin.Rda")

## count non-missing obs ## 

sum(!is.na(af.muni.fin$mgr_npov_2022))
sum(!is.na(af.muni.fin$mgr_hpov_2022))
sum(!is.na(af.muni.fin$mgr_hpov_cc_2022))

sum(!is.na(af.muni.fin$rb_npov_2022))
sum(!is.na(af.muni.fin$rb_hpov_2022))
sum(!is.na(af.muni.fin$rb_hpov_cc_2022))

sum(!is.na(af.muni.fin$mpv_npov_2022))
sum(!is.na(af.muni.fin$mpv_hpov_2022))
sum(!is.na(af.muni.fin$mpv_hpov_cc_2022))

## function to extract the results ##

get_res <- function(i,j){
  
  inres <- qread(paste0(data_path,
                        "202_muni_results_supp_",
                        i,
                        ".qs"))
  
  inres.fin <- inres[[j]]
  return(inres.fin)
  
}

muni.noimp.res1 <- list()
muni.noimp.res2 <- list()
muni.noimp.res3 <- list()
muni.noimp.res4 <- list()
muni.noimp.res5 <- list()
muni.noimp.res6 <- list()
muni.noimp.res7 <- list()
muni.noimp.res8 <- list()
muni.noimp.res9 <- list()

muni.noimp.res1 <- mapply(get_res,
                          1:100,
                          10,
                          SIMPLIFY = F)

save(muni.noimp.res1, 
     file = "203f_m1_muni_noimp.Rda")

muni.noimp.res2 <- mapply(get_res,
                          1:100,
                          11,
                          SIMPLIFY = F)

save(muni.noimp.res2, 
     file = "203f_m2_muni_noimp.Rda")

muni.noimp.res3 <- mapply(get_res,
                          1:100,
                          12,
                          SIMPLIFY = F)

save(muni.noimp.res3, 
     file = "203f_m3_muni_noimp.Rda")

muni.noimp.res4 <- mapply(get_res,
                          1:100,
                          13,
                          SIMPLIFY = F)

save(muni.noimp.res4, 
     file = "203f_m4_muni_noimp.Rda")

muni.noimp.res5 <- mapply(get_res,
                          1:100,
                          14,
                          SIMPLIFY = F)

save(muni.noimp.res5, 
     file = "203f_m5_muni_noimp.Rda")

muni.noimp.res6 <- mapply(get_res,
                          1:100,
                          15,
                          SIMPLIFY = F)

save(muni.noimp.res6, 
     file = "203f_m6_muni_noimp.Rda")

muni.noimp.res7 <- mapply(get_res,
                          1:100,
                          16,
                          SIMPLIFY = F)

save(muni.noimp.res7, 
     file = "203f_m7_muni_noimp.Rda")

muni.noimp.res8 <- mapply(get_res,
                          1:100,
                          17,
                          SIMPLIFY = F)

save(muni.noimp.res8, 
     file = "203f_m8_muni_noimp.Rda")

muni.noimp.res9 <- mapply(get_res,
                          1:100,
                          18,
                          SIMPLIFY = F)

save(muni.noimp.res9, 
     file = "203f_m9_muni_noimp.Rda")

## compile the results ## 

ex.msm.res <- function(i, inlist) {
  
  #if(is.na(inlist[[i]][1])){
  #  return(NA)
  #}

  res <- as.data.frame(inlist[[i]][1])
  
  return(res)
  
}

ex.noadj.res <- function(i, inlist) {
  
  #if(is.na(inlist[[i]][3])){
  #  return(NA)
  #}
  
  res <- as.data.frame(inlist[[i]][3])
  
  return(res)
  
}

ex.adl.res <- function(i, inlist) {
  
  #if(is.na(inlist[[i]][4])){
  #  return(NA)
  #}
  
  res <- as.data.frame(inlist[[i]][4])
  
  return(res)
  
}

ex.fe.res <- function(i, inlist) {
  
  #if(is.na(inlist[[i]][5])){
  #  return(NA)
  #}
  
  res <- as.data.frame(inlist[[i]][5])
  
  return(res)
  
}

## load the results ##

load("203f_m1_muni_noimp.Rda")
load("203f_m2_muni_noimp.Rda")
load("203f_m3_muni_noimp.Rda")
load("203f_m4_muni_noimp.Rda")
load("203f_m5_muni_noimp.Rda")
load("203f_m6_muni_noimp.Rda")
load("203f_m7_muni_noimp.Rda")
load("203f_m8_muni_noimp.Rda")
load("203f_m9_muni_noimp.Rda")

##
## msm results ## 
##

m1.msm.noimp.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.noimp.res1)

m1.msm.noimp.res.ov <- bind_rows(m1.msm.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist")) 

m2.msm.noimp.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.noimp.res2)

m2.msm.noimp.res.ov <- bind_rows(m2.msm.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist")) 

m3.msm.noimp.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.noimp.res3)

m3.msm.noimp.res.ov <- bind_rows(m3.msm.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist")) 

m4.msm.noimp.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.noimp.res4)

m4.msm.noimp.res.ov <- bind_rows(m4.msm.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist"))

m5.msm.noimp.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.noimp.res5)

m5.msm.noimp.res.ov <- bind_rows(m5.msm.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist"))

m6.msm.noimp.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.noimp.res6)

m6.msm.noimp.res.ov <- bind_rows(m6.msm.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist"))

m7.msm.noimp.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.noimp.res7)

m7.msm.noimp.res.ov <- bind_rows(m7.msm.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist"))

m8.msm.noimp.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.noimp.res8)

m8.msm.noimp.res.ov <- bind_rows(m8.msm.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist"))

m9.msm.noimp.res <- lapply(seq(1,100,by=1),
                     ex.msm.res,
                     inlist = muni.noimp.res9)

m9.msm.noimp.res.ov <- bind_rows(m9.msm.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist"))


##
## noadj results ## 
##

m1.noadj.noimp.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.noimp.res1)

m1.noadj.noimp.res.ov <- bind_rows(m1.noadj.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist"))

m2.noadj.noimp.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.noimp.res2)

m2.noadj.noimp.res.ov <- bind_rows(m2.noadj.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist"))

m3.noadj.noimp.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.noimp.res3)

m3.noadj.noimp.res.ov <- bind_rows(m3.noadj.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist"))

m4.noadj.noimp.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.noimp.res4)

m4.noadj.noimp.res.ov <- bind_rows(m4.noadj.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist"))

m5.noadj.noimp.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.noimp.res5)

m5.noadj.noimp.res.ov <- bind_rows(m5.noadj.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist"))

m6.noadj.noimp.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.noimp.res6)

m6.noadj.noimp.res.ov <- bind_rows(m6.noadj.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist"))

m7.noadj.noimp.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.noimp.res7)

m7.noadj.noimp.res.ov <- bind_rows(m7.noadj.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist"))

m8.noadj.noimp.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.noimp.res8)

m8.noadj.noimp.res.ov <- bind_rows(m8.noadj.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist"))

m9.noadj.noimp.res <- lapply(seq(1,100,by=1),
                       ex.noadj.res,
                       inlist = muni.noimp.res9)

m9.noadj.noimp.res.ov <- bind_rows(m9.noadj.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "thist"))


##
## adl results ## 
##

m1.adl.noimp.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.noimp.res1)

m1.adl.noimp.res.ov <- bind_rows(m1.adl.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_2022_fin_st"))

m2.adl.noimp.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.noimp.res2)

m2.adl.noimp.res.ov <- bind_rows(m2.adl.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_2022_fin_st"))

m3.adl.noimp.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.noimp.res3)

m3.adl.noimp.res.ov <- bind_rows(m3.adl.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_2022_fin_st"))

m4.adl.noimp.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.noimp.res4)

m4.adl.noimp.res.ov <- bind_rows(m4.adl.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_2022_fin_st"))

m5.adl.noimp.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.noimp.res5)

m5.adl.noimp.res.ov <- bind_rows(m5.adl.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_2022_fin_st"))

m6.adl.noimp.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.noimp.res6)

m6.adl.noimp.res.ov <- bind_rows(m6.adl.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_2022_fin_st"))

m7.adl.noimp.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.noimp.res7)

m7.adl.noimp.res.ov <- bind_rows(m7.adl.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_2022_fin_st"))

m8.adl.noimp.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.noimp.res8)

m8.adl.noimp.res.ov <- bind_rows(m8.adl.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_2022_fin_st"))

m9.adl.noimp.res <- lapply(seq(1,100,by=1),
                     ex.adl.res,
                     inlist = muni.noimp.res9)

m9.adl.noimp.res.ov <- bind_rows(m9.adl.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_2022_fin_st"))

## fixed effects ##

m1.fe.noimp.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.noimp.res1)

m1.fe.noimp.res.ov <- bind_rows(m1.fe.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_fin"))

m2.fe.noimp.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.noimp.res2)

m2.fe.noimp.res.ov <- bind_rows(m2.fe.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_fin"))

m3.fe.noimp.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.noimp.res3)

m3.fe.noimp.res.ov <- bind_rows(m3.fe.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_fin"))

m4.fe.noimp.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.noimp.res4)

m4.fe.noimp.res.ov <- bind_rows(m4.fe.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_fin"))

m5.fe.noimp.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.noimp.res5)

m5.fe.noimp.res.ov <- bind_rows(m5.fe.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_fin"))

m6.fe.noimp.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.noimp.res6)

m6.fe.noimp.res.ov <- bind_rows(m6.fe.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_fin"))

m7.fe.noimp.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.noimp.res7)

m7.fe.noimp.res.ov <- bind_rows(m7.fe.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_fin"))

m8.fe.noimp.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.noimp.res8)

m8.fe.noimp.res.ov <- bind_rows(m8.fe.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_fin"))

m9.fe.noimp.res <- lapply(seq(1,100,by=1),
                    ex.fe.res,
                    inlist = muni.noimp.res9)

m9.fe.noimp.res.ov <- bind_rows(m9.fe.noimp.res) %>%
  rownames_to_column(var = "term") %>%
  filter(str_starts(term, "zri_fin"))

## plot the results ##

msm.noimp.res.in <- data.frame(var = c("Median gross rent (nat)",
                                 "Median gross rent (hpov)",
                                 "Median gross rent (hpov cc)",
                                 "Rent burden (nat)",
                                 "Rent burden (hpov)",
                                 "Rent burden (hpov cc)",
                                 "MPV (nat)",
                                 "MPV (hpov)",
                                 "MPV (hpov cc)"),
                         coeff = c(mean(m1.msm.noimp.res.ov$Estimate, na.rm=T),
                                   mean(m2.msm.noimp.res.ov$Estimate, na.rm=T),
                                   mean(m3.msm.noimp.res.ov$Estimate, na.rm=T),
                                   mean(m4.msm.noimp.res.ov$Estimate, na.rm=T),
                                   mean(m5.msm.noimp.res.ov$Estimate, na.rm=T),
                                   mean(m6.msm.noimp.res.ov$Estimate, na.rm=T),
                                   mean(m7.msm.noimp.res.ov$Estimate, na.rm=T),
                                   mean(m8.msm.noimp.res.ov$Estimate, na.rm=T),
                                   mean(m9.msm.noimp.res.ov$Estimate, na.rm=T)),
                         lb = c(quantile(m1.msm.noimp.res.ov$Estimate,0.025, na.rm=T),
                                quantile(m2.msm.noimp.res.ov$Estimate,0.025, na.rm=T),
                                quantile(m3.msm.noimp.res.ov$Estimate,0.025, na.rm=T),
                                quantile(m4.msm.noimp.res.ov$Estimate,0.025, na.rm=T),
                                quantile(m5.msm.noimp.res.ov$Estimate,0.025, na.rm=T),
                                quantile(m6.msm.noimp.res.ov$Estimate,0.025, na.rm=T),
                                quantile(m7.msm.noimp.res.ov$Estimate,0.025, na.rm=T),
                                quantile(m8.msm.noimp.res.ov$Estimate,0.025, na.rm=T),
                                quantile(m9.msm.noimp.res.ov$Estimate,0.025, na.rm=T)),
                         ub = c(quantile(m1.msm.noimp.res.ov$Estimate,0.975, na.rm=T),
                                quantile(m2.msm.noimp.res.ov$Estimate,0.975, na.rm=T),
                                quantile(m3.msm.noimp.res.ov$Estimate,0.975, na.rm=T),
                                quantile(m4.msm.noimp.res.ov$Estimate,0.975, na.rm=T),
                                quantile(m5.msm.noimp.res.ov$Estimate,0.975, na.rm=T),
                                quantile(m6.msm.noimp.res.ov$Estimate,0.975, na.rm=T),
                                quantile(m7.msm.noimp.res.ov$Estimate,0.975, na.rm=T),
                                quantile(m8.msm.noimp.res.ov$Estimate,0.975, na.rm=T),
                                quantile(m9.msm.noimp.res.ov$Estimate,0.975, na.rm=T)))

## no adj ## 

noadj.noimp.res.in <- data.frame(var = c("Median gross rent (nat)",
                                   "Median gross rent (hpov)",
                                   "Median gross rent (hpov cc)",
                                   "Rent burden (nat)",
                                   "Rent burden (hpov)",
                                   "Rent burden (hpov cc)",
                                   "MPV (nat)",
                                   "MPV (hpov)",
                                   "MPV (hpov cc)"),
                           coeff = c(mean(m1.noadj.noimp.res.ov$Estimate, na.rm=T),
                                     mean(m2.noadj.noimp.res.ov$Estimate, na.rm=T),
                                     mean(m3.noadj.noimp.res.ov$Estimate, na.rm=T),
                                     mean(m4.noadj.noimp.res.ov$Estimate, na.rm=T),
                                     mean(m5.noadj.noimp.res.ov$Estimate, na.rm=T),
                                     mean(m6.noadj.noimp.res.ov$Estimate, na.rm=T),
                                     mean(m7.noadj.noimp.res.ov$Estimate, na.rm=T),
                                     mean(m8.noadj.noimp.res.ov$Estimate, na.rm=T),
                                     mean(m9.noadj.noimp.res.ov$Estimate, na.rm=T)),
                           lb = c(quantile(m1.noadj.noimp.res.ov$Estimate,0.025, na.rm=T),
                                  quantile(m2.noadj.noimp.res.ov$Estimate,0.025, na.rm=T),
                                  quantile(m3.noadj.noimp.res.ov$Estimate,0.025, na.rm=T),
                                  quantile(m4.noadj.noimp.res.ov$Estimate,0.025, na.rm=T),
                                  quantile(m5.noadj.noimp.res.ov$Estimate,0.025, na.rm=T),
                                  quantile(m6.noadj.noimp.res.ov$Estimate,0.025, na.rm=T),
                                  quantile(m7.noadj.noimp.res.ov$Estimate,0.025, na.rm=T),
                                  quantile(m8.noadj.noimp.res.ov$Estimate,0.025, na.rm=T),
                                  quantile(m9.noadj.noimp.res.ov$Estimate,0.025, na.rm=T)),
                           ub = c(quantile(m1.noadj.noimp.res.ov$Estimate,0.975, na.rm=T),
                                  quantile(m2.noadj.noimp.res.ov$Estimate,0.975, na.rm=T),
                                  quantile(m3.noadj.noimp.res.ov$Estimate,0.975, na.rm=T),
                                  quantile(m4.noadj.noimp.res.ov$Estimate,0.975, na.rm=T),
                                  quantile(m5.noadj.noimp.res.ov$Estimate,0.975, na.rm=T),
                                  quantile(m6.noadj.noimp.res.ov$Estimate,0.975, na.rm=T),
                                  quantile(m7.noadj.noimp.res.ov$Estimate,0.975, na.rm=T),
                                  quantile(m8.noadj.noimp.res.ov$Estimate,0.975, na.rm=T),
                                  quantile(m9.noadj.noimp.res.ov$Estimate,0.975, na.rm=T)))

## adl ## 

adl.noimp.res.in <- data.frame(var = c("Median gross rent (nat)",
                                 "Median gross rent (hpov)",
                                 "Median gross rent (hpov cc)",
                                 "Rent burden (nat)",
                                 "Rent burden (hpov)",
                                 "Rent burden (hpov cc)",
                                 "MPV (nat)",
                                 "MPV (hpov)",
                                 "MPV (hpov cc)"),
                         coeff = c(mean(m1.adl.noimp.res.ov$Estimate, na.rm=T),
                                   mean(m2.adl.noimp.res.ov$Estimate, na.rm=T),
                                   mean(m3.adl.noimp.res.ov$Estimate, na.rm=T),
                                   mean(m4.adl.noimp.res.ov$Estimate, na.rm=T),
                                   mean(m5.adl.noimp.res.ov$Estimate, na.rm=T),
                                   mean(m6.adl.noimp.res.ov$Estimate, na.rm=T),
                                   mean(m7.adl.noimp.res.ov$Estimate, na.rm=T),
                                   mean(m8.adl.noimp.res.ov$Estimate, na.rm=T),
                                   mean(m9.adl.noimp.res.ov$Estimate, na.rm=T)),
                         lb = c(quantile(m1.adl.noimp.res.ov$Estimate,0.025, na.rm=T),
                                quantile(m2.adl.noimp.res.ov$Estimate,0.025, na.rm=T),
                                quantile(m3.adl.noimp.res.ov$Estimate,0.025, na.rm=T),
                                quantile(m4.adl.noimp.res.ov$Estimate,0.025, na.rm=T),
                                quantile(m5.adl.noimp.res.ov$Estimate,0.025, na.rm=T),
                                quantile(m6.adl.noimp.res.ov$Estimate,0.025, na.rm=T),
                                quantile(m7.adl.noimp.res.ov$Estimate,0.025, na.rm=T),
                                quantile(m8.adl.noimp.res.ov$Estimate,0.025, na.rm=T),
                                quantile(m9.adl.noimp.res.ov$Estimate,0.025, na.rm=T)),
                         ub = c(quantile(m1.adl.noimp.res.ov$Estimate,0.975, na.rm=T),
                                quantile(m2.adl.noimp.res.ov$Estimate,0.975, na.rm=T),
                                quantile(m3.adl.noimp.res.ov$Estimate,0.975, na.rm=T),
                                quantile(m4.adl.noimp.res.ov$Estimate,0.975, na.rm=T),
                                quantile(m5.adl.noimp.res.ov$Estimate,0.975, na.rm=T),
                                quantile(m6.adl.noimp.res.ov$Estimate,0.975, na.rm=T),
                                quantile(m7.adl.noimp.res.ov$Estimate,0.975, na.rm=T),
                                quantile(m8.adl.noimp.res.ov$Estimate,0.975, na.rm=T),
                                quantile(m9.adl.noimp.res.ov$Estimate,0.975, na.rm=T)))


fe.noimp.res.in <- data.frame(var = c("Median gross rent (nat)",
                                "Median gross rent (hpov)",
                                "Median gross rent (hpov cc)",
                                "Rent burden (nat)",
                                "Rent burden (hpov)",
                                "Rent burden (hpov cc)",
                                "MPV (nat)",
                                "MPV (hpov)",
                                "MPV (hpov cc)"),
                        coeff = c(mean(m1.fe.noimp.res.ov$Estimate, na.rm=T),
                                  mean(m2.fe.noimp.res.ov$Estimate, na.rm=T),
                                  mean(m3.fe.noimp.res.ov$Estimate, na.rm=T),
                                  mean(m4.fe.noimp.res.ov$Estimate, na.rm=T),
                                  mean(m5.fe.noimp.res.ov$Estimate, na.rm=T),
                                  mean(m6.fe.noimp.res.ov$Estimate, na.rm=T),
                                  mean(m7.fe.noimp.res.ov$Estimate, na.rm=T),
                                  mean(m8.fe.noimp.res.ov$Estimate, na.rm=T),
                                  mean(m9.fe.noimp.res.ov$Estimate, na.rm=T)),
                        lb = c(quantile(m1.fe.noimp.res.ov$Estimate,0.025, na.rm=T),
                               quantile(m2.fe.noimp.res.ov$Estimate,0.025, na.rm=T),
                               quantile(m3.fe.noimp.res.ov$Estimate,0.025, na.rm=T),
                               quantile(m4.fe.noimp.res.ov$Estimate,0.025, na.rm=T),
                               quantile(m5.fe.noimp.res.ov$Estimate,0.025, na.rm=T),
                               quantile(m6.fe.noimp.res.ov$Estimate,0.025, na.rm=T),
                               quantile(m7.fe.noimp.res.ov$Estimate,0.025, na.rm=T),
                               quantile(m8.fe.noimp.res.ov$Estimate,0.025, na.rm=T),
                               quantile(m9.fe.noimp.res.ov$Estimate,0.025, na.rm=T)),
                        ub = c(quantile(m1.fe.noimp.res.ov$Estimate,0.975, na.rm=T),
                               quantile(m2.fe.noimp.res.ov$Estimate,0.975, na.rm=T),
                               quantile(m3.fe.noimp.res.ov$Estimate,0.975, na.rm=T),
                               quantile(m4.fe.noimp.res.ov$Estimate,0.975, na.rm=T),
                               quantile(m5.fe.noimp.res.ov$Estimate,0.975, na.rm=T),
                               quantile(m6.fe.noimp.res.ov$Estimate,0.975, na.rm=T),
                               quantile(m7.fe.noimp.res.ov$Estimate,0.975, na.rm=T),
                               quantile(m8.fe.noimp.res.ov$Estimate,0.975, na.rm=T),
                               quantile(m9.fe.noimp.res.ov$Estimate,0.975, na.rm=T)))


## set 1: median rents ## 

msm.fin.noimp.res1 <- msm.noimp.res.in %>%
  filter(var %in% c("Median gross rent (nat)",
                    "Median gross rent (hpov)",
                    "Median gross rent (hpov cc)")) %>%
  mutate(group = case_when(var == "Median gross rent (nat)" ~ "Above-average\n(N=1,082)",
                           var == "Median gross rent (hpov)" ~ "High-poverty\n(N=233)",
                           var == "Median gross rent (hpov cc)" ~ "High-poverty\ncentral city\n(N=189)"),
         var = case_when(var == "Median gross rent (nat)" ~ "Above-average 1 (N=1,082)",
                         var == "Median gross rent (hpov)" ~ "High-poverty 1 (N=233)",
                         var == "Median gross rent (hpov cc)" ~ "High-poverty cc 1 (N=189)"),
         lb = lb,
         ub = ub,
         Model = "MSM")

noadj.fin.noimp.res1 <- noadj.noimp.res.in %>%
  filter(var %in% c("Median gross rent (nat)",
                    "Median gross rent (hpov)",
                    "Median gross rent (hpov cc)")) %>%
  mutate(group = case_when(var == "Median gross rent (nat)" ~ "Above-average\n(N=1,082)",
                           var == "Median gross rent (hpov)" ~ "High-poverty\n(N=233)",
                           var == "Median gross rent (hpov cc)" ~ "High-poverty\ncentral city\n(N=189)"),
         var = case_when(var == "Median gross rent (nat)" ~ "Above-average 2 (N=1,082)",
                         var == "Median gross rent (hpov)" ~ "High-poverty 2 (N=233)",
                         var == "Median gross rent (hpov cc)" ~ "High-poverty cc 2 (N=189)"),
         lb = lb,
         ub = ub,
         Model = "NoAdj")

adl.fin.noimp.res1 <- adl.noimp.res.in %>%
  filter(var %in% c("Median gross rent (nat)",
                    "Median gross rent (hpov)",
                    "Median gross rent (hpov cc)")) %>%
  mutate(group = case_when(var == "Median gross rent (nat)" ~ "Above-average\n(N=1,082)",
                           var == "Median gross rent (hpov)" ~ "High-poverty\n(N=233)",
                           var == "Median gross rent (hpov cc)" ~ "High-poverty\ncentral city\n(N=189)"),
         var = case_when(var == "Median gross rent (nat)" ~ "Above-average 3 (N=1,082)",
                         var == "Median gross rent (hpov)" ~ "High-poverty 3 (N=233)",
                         var == "Median gross rent (hpov cc)" ~ "High-poverty cc 3 (N=189)"),
         lb = lb,
         ub = ub,
         Model = "ADL")

fe.fin.noimp.res1 <- fe.noimp.res.in %>%
  filter(var %in% c("Median gross rent (nat)",
                    "Median gross rent (hpov)",
                    "Median gross rent (hpov cc)")) %>%
  mutate(group = case_when(var == "Median gross rent (nat)" ~ "Above-average\n(N=1,082)",
                           var == "Median gross rent (hpov)" ~ "High-poverty\n(N=233)",
                           var == "Median gross rent (hpov cc)" ~ "High-poverty\ncentral city\n(N=189)"),
         var = case_when(var == "Median gross rent (nat)" ~ "Above-average 4 (N=1,082)",
                         var == "Median gross rent (hpov)" ~ "High-poverty 4 (N=233)",
                         var == "Median gross rent (hpov cc)" ~ "High-poverty cc 4 (N=189)"),
         lb = lb,
         ub = ub,
         Model = "FE")

fin.noimp.res1 <- rbind(msm.fin.noimp.res1,
                  noadj.fin.noimp.res1,
                  adl.fin.noimp.res1,
                  fe.fin.noimp.res1)

fin.noimp.res1$Model <- factor(fin.noimp.res1$Model,
                         levels = c("MSM",
                                    "NoAdj",
                                    "ADL",
                                    "FE"))


ggplot(fin.noimp.res1, aes(x=var, y=coeff, color=Model)) + 
  geom_pointrange(aes(ymin=lb, ymax=ub)) + 
  geom_hline(yintercept = 0) +
  theme_bw(base_size = 14) + 
  theme(
    axis.text.x = element_blank(), # Hide individual x-axis labels
    axis.title.x = element_blank(),
    axis.ticks.x = element_blank(),
    panel.spacing = unit(0, "lines"), 
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text.x = element_text(size = 14, face = "bold"), # Adjust group label text
    legend.title = element_text(size = 14, face = "bold"),
    legend.text = element_text(size = 14, face = "bold"),
  ) + 
  facet_wrap(~group, strip.position = "bottom", scales = "free_x") +
  ggtitle("") + 
  xlab("") + 
  ylab("Estimate ($)") + 
  scale_colour_manual(values = c("royalblue","forestgreen","red","orange"))


## set 2: rent burden ##

msm.fin.noimp.res2 <- msm.noimp.res.in %>%
  filter(var %in% c("Rent burden (nat)",
                    "Rent burden (hpov)",
                    "Rent burden (hpov cc)")) %>%
  mutate(group = case_when(var == "Rent burden (nat)" ~ "Above-average\n(N=1,082)",
                           var == "Rent burden (hpov)" ~ "High-poverty\n(N=234)",
                           var == "Rent burden (hpov cc)" ~ "High-poverty\ncentral city\n(N=189)"),
         var = case_when(var == "Rent burden (nat)" ~ "Above-average 1 (N=1,082)",
                         var == "Rent burden (hpov)" ~ "High-poverty 1 (N=234)",
                         var == "Rent burden (hpov cc)" ~ "High-poverty cc 1 (N=189)"),
         lb = lb,
         ub = ub,
         Model = "MSM")

noadj.fin.noimp.res2 <- noadj.noimp.res.in %>%
  filter(var %in% c("Rent burden (nat)",
                    "Rent burden (hpov)",
                    "Rent burden (hpov cc)")) %>%
  mutate(group = case_when(var == "Rent burden (nat)" ~ "Above-average\n(N=1,082)",
                           var == "Rent burden (hpov)" ~ "High-poverty\n(N=234)",
                           var == "Rent burden (hpov cc)" ~ "High-poverty\ncentral city\n(N=189)"),
         var = case_when(var == "Rent burden (nat)" ~ "Above-average 2 (N=1,082)",
                         var == "Rent burden (hpov)" ~ "High-poverty 2 (N=234)",
                         var == "Rent burden (hpov cc)" ~ "High-poverty cc 2 (N=189)"),
         lb = lb,
         ub = ub,
         Model = "NoAdj")

adl.fin.noimp.res2 <- adl.noimp.res.in %>%
  filter(var %in% c("Rent burden (nat)",
                    "Rent burden (hpov)",
                    "Rent burden (hpov cc)")) %>%
  mutate(group = case_when(var == "Rent burden (nat)" ~ "Above-average\n(N=1,082)",
                           var == "Rent burden (hpov)" ~ "High-poverty\n(N=234)",
                           var == "Rent burden (hpov cc)" ~ "High-poverty\ncentral city\n(N=189)"),
         var = case_when(var == "Rent burden (nat)" ~ "Above-average 3 (N=1,082)",
                         var == "Rent burden (hpov)" ~ "High-poverty 3 (N=234)",
                         var == "Rent burden (hpov cc)" ~ "High-poverty cc 3 (N=189)"),
         lb = lb,
         ub = ub,
         Model = "ADL")

fe.fin.noimp.res2 <- fe.noimp.res.in %>%
  filter(var %in% c("Rent burden (nat)",
                    "Rent burden (hpov)",
                    "Rent burden (hpov cc)")) %>%
  mutate(group = case_when(var == "Rent burden (nat)" ~ "Above-average\n(N=1,082)",
                           var == "Rent burden (hpov)" ~ "High-poverty\n(N=234)",
                           var == "Rent burden (hpov cc)" ~ "High-poverty\ncentral city\n(N=189)"),
         var = case_when(var == "Rent burden (nat)" ~ "Above-average 4 (N=1,082)",
                         var == "Rent burden (hpov)" ~ "High-poverty 4 (N=234)",
                         var == "Rent burden (hpov cc)" ~ "High-poverty cc 4 (N=189)"),
         lb = lb,
         ub = ub,
         Model = "FE")

fin.noimp.res2 <- rbind(msm.fin.noimp.res2,
                  noadj.fin.noimp.res2,
                  adl.fin.noimp.res2,
                  fe.fin.noimp.res2)

fin.noimp.res2$Model <- factor(fin.noimp.res2$Model,
                         levels = c("MSM",
                                    "NoAdj",
                                    "ADL",
                                    "FE"))


ggplot(fin.noimp.res2, aes(x=var, y=coeff, color=Model)) + 
  geom_pointrange(aes(ymin=lb, ymax=ub)) + 
  geom_hline(yintercept = 0) +
  theme_bw(base_size = 14) + 
  theme(
    axis.text.x = element_blank(), # Hide individual x-axis labels
    axis.title.x = element_blank(),
    axis.ticks.x = element_blank(),
    panel.spacing = unit(0, "lines"), 
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text.x = element_text(size = 14, face = "bold"), # Adjust group label text
    legend.title = element_text(size = 14, face = "bold"),
    legend.text = element_text(size = 14, face = "bold"),
  ) + 
  facet_wrap(~group, strip.position = "bottom", scales = "free_x") +
  ggtitle("") + 
  xlab("") + 
  ylab("Estimate") + 
  scale_colour_manual(values = c("royalblue","forestgreen","red","orange"))

## set 3: MPV in impoverished neighborhoods ##

msm.fin.noimp.res3 <- msm.noimp.res.in %>%
  filter(var %in% c("MPV (nat)",
                    "MPV (hpov)",
                    "MPV (hpov cc)")) %>%
  mutate(group = case_when(var == "MPV (nat)" ~ "Above-average\n(N=1,081)",
                           var == "MPV (hpov)" ~ "High-poverty\n(N=227)",
                           var == "MPV (hpov cc)" ~ "High-poverty\ncentral city\n(N=184)"),
         var = case_when(var == "MPV (nat)" ~ "Above-average 1 (N=1,081)",
                         var == "MPV (hpov)" ~ "High-poverty 1 (N=227)",
                         var == "MPV (hpov cc)" ~ "High-poverty cc 1 (N=184)"),
         lb = lb,
         ub = ub,
         Model = "MSM")

noadj.fin.noimp.res3 <- noadj.noimp.res.in %>%
  filter(var %in% c("MPV (nat)",
                    "MPV (hpov)",
                    "MPV (hpov cc)")) %>%
  mutate(group = case_when(var == "MPV (nat)" ~ "Above-average\n(N=1,081)",
                           var == "MPV (hpov)" ~ "High-poverty\n(N=227)",
                           var == "MPV (hpov cc)" ~ "High-poverty\ncentral city\n(N=184)"),
         var = case_when(var == "MPV (nat)" ~ "Above-average 2 (N=1,081)",
                         var == "MPV (hpov)" ~ "High-poverty 2 (N=227)",
                         var == "MPV (hpov cc)" ~ "High-poverty cc 2 (N=184)"),
         lb = lb,
         ub = ub,
         Model = "NoAdj")

adl.fin.noimp.res3 <- adl.noimp.res.in %>%
  filter(var %in% c("MPV (nat)",
                    "MPV (hpov)",
                    "MPV (hpov cc)")) %>%
  mutate(group = case_when(var == "MPV (nat)" ~ "Above-average\n(N=1,081)",
                           var == "MPV (hpov)" ~ "High-poverty\n(N=227)",
                           var == "MPV (hpov cc)" ~ "High-poverty\ncentral city\n(N=184)"),
         var = case_when(var == "MPV (nat)" ~ "Above-average 3 (N=1,081)",
                         var == "MPV (hpov)" ~ "High-poverty 3 (N=227)",
                         var == "MPV (hpov cc)" ~ "High-poverty cc 3 (N=184)"),
         lb = lb,
         ub = ub,
         Model = "ADL")


fe.fin.noimp.res3 <- fe.noimp.res.in %>%
  filter(var %in% c("MPV (nat)",
                    "MPV (hpov)",
                    "MPV (hpov cc)")) %>%
  mutate(group = case_when(var == "MPV (nat)" ~ "Above-average\n(N=1,081)",
                           var == "MPV (hpov)" ~ "High-poverty\n(N=227)",
                           var == "MPV (hpov cc)" ~ "High-poverty\ncentral city\n(N=184)"),
         var = case_when(var == "MPV (nat)" ~ "Above-average 4 (N=1,081)",
                         var == "MPV (hpov)" ~ "High-poverty 4 (N=227)",
                         var == "MPV (hpov cc)" ~ "High-poverty cc 4 (N=184)"),
         lb = lb,
         ub = ub,
         Model = "FE")

fin.noimp.res3 <- rbind(msm.fin.noimp.res3,
                  noadj.fin.noimp.res3,
                  adl.fin.noimp.res3,
                  fe.fin.noimp.res3)

fin.noimp.res3$Model <- factor(fin.noimp.res3$Model,
                         levels = c("MSM",
                                    "NoAdj",
                                    "ADL",
                                    "FE"))

ggplot(fin.noimp.res3, aes(x=var, y=coeff, color=Model)) + 
  geom_pointrange(aes(ymin=lb, ymax=ub)) + 
  geom_hline(yintercept = 0) +
  theme_bw(base_size = 14) + 
  theme(
    axis.text.x = element_blank(), # Hide individual x-axis labels
    axis.title.x = element_blank(),
    axis.ticks.x = element_blank(),
    panel.spacing = unit(0, "lines"), 
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text.x = element_text(size = 14, face = "bold"), # Adjust group label text
    legend.title = element_text(size = 14, face = "bold"),
    legend.text = element_text(size = 14, face = "bold"),
  ) + 
  facet_wrap(~group, strip.position = "bottom", scales = "free_x") +
  ggtitle("") + 
  xlab("") + 
  ylab("Estimate ($)") + 
  scale_colour_manual(values = c("royalblue","forestgreen","red","orange"))


## END OF PROGRAM ##

#sink()
